tic
% Построение реализаций при шумовом воздействии
% 	Фролов Д.А. ЭР-23-09
clc       % очистка экрана
close all % закрытие всех лишних окон
L=100e-6; % Катушка индуктивности, Гн
global L  % Глобальная переменная, чтобы использовать в функции f
R=47;	  % Сопротивление резистора в цепи, Ом
global R  % Глобальная переменная, чтобы использовать в функции f
A=1;      % амплитуда входного гармонического воздействия, В
global A  % Глобальная переменная, чтобы использовать в функции f

Nt=100; 		        % количество отсчетов по времени
tmin=0;         		% начало временного интервала анализа
tmax=10e-6;       		% конец временного анализа, 75мкс, если взять 50, то на малых частотах получим скачок
dt=(tmax-tmin)/Nt; 		% шаг по времени, 7.5 нс
t=(0:Nt-1)*dt;  		% задаю вектор времени
global dt 		        % Глобальная переменная, чтобы использовать в функции f

x0=0;             		% начальное значение

sko=1;             		% СКО шума
m=0;               		% математическое ожидание шума
noise=randn(1,Nt)*sko+m;% заполняем вектор-реализацию шума
global noise  			% Глобальная переменная, чтобы использовать в функции f

% Для решения ДУ в мат. пакете Octave есть функция lsode(f,x0,t)
% f  - имя функции, которая описывает ДУ, имеет два входных параметра
% х0 - начальное условие в первый момент времени
% t  - вектор, который содержит отсчеты времени. Первый элемент - для начального условия
function xdot=f(x,t)
	global L     % используем заранее объявленную индуктивность
	global R     % используем заранее объявленное сопротивление
	global A     % используем заранее объявленную амплитуду
	global noise % используем заранее объявленный шум
	global dt    % используем заранее объявленный шаг по частоте
	xdot=-2*R/L*x+noise(round(t/dt)+1)/L;
endfunction

cur=lsode('f',x0,t);              % рещаем ДУ, находим реализацию на выходе
plot(t/1e-6,noise,t/1e-6,cur*R)   % строим реализацию входного и выходного шума
xlabel('t, us')          		  % подпись оси ОХ
ylabel('Uin, Uout, V')            % подпись оси OY
title('OSC AWGN')     		      % Название графика
grid on                           % вывести сетку
box off                           % оставить только оси
legend('Uin','Uout')              % вывести легенду
ylim([min(noise),max(noise)])     % границы отображения по оси OY

sko=0.05;    
noise=randn(1,Nt)*sko+m;
for i=1:Nt
	noise(i)=1*sin(2e5*2*pi*t(i))+noise(i);
	end
figure(2)
cur=lsode('f',x0,t);              % рещаем ДУ, находим реализацию на выходе
plot(t/1e-6,noise,t/1e-6,cur*R)   % строим реализацию входного и выходного шума
xlabel('t, us')          		  % подпись оси ОХ
								  % подпись оси OY
ylabel('Uin with AWGN, Uout without AWGN, V')            
title('OSC sin and AWGN')     	  % Название графика
grid on                           % вывести сетку
box off                           % оставить только оси
legend('Uin','Uout')              % вывести легенду
ylim([min(noise),max(noise)])     % границы отображения по оси OY
toc